classdef GetXe131T1Rho < Table1.NMRGJobs.NMRGJob
    methods
        function obj = GetXe131T1Rho(nmrg)
            obj@Table1.NMRGJobs.NMRGJob(nmrg, 'GetXe131T1Rho');
            obj.defineParameter();
        end
        
        function obj = defineParameter(obj)
            obj.paramParser=inputParser;
            obj.paramParser.addParameter('duration', 60, @(x) x>=0);
            obj.paramParser.addParameter('durationDecay', 120, @(x) x>=0);
            obj.paramParser.addParameter('amplitude', 1.0, @(x) x>0.0 && x<=2.0);
            obj.paramParser.addParameter('detuning', 0.0, @isnumeric);
            obj.paramParser.addParameter('preDrivingTime', 10.0, @(x) x>0);
            obj.paramParser.addParameter('calibrationField', 0.02, @(x) x>0.0);
            obj.paramParser.addParameter('sign', 1, @(x) ismember(x, [-1 1]));
            obj.paramParser.addParameter('rate', 50, @(x) x>0);
            obj.paramParser.addParameter('isPlot', true, @islogical);
            obj.paramParser.addParameter('description', 'No Description.', @ischar);
        end
        %%
        function obj = preprocess(obj)
            nmrg = obj.boss;

            nmrg.UnlockXe129(); nmrg.UnlockXe131();
            nmrg.setRFMagnetometer('frequency', nmrg.SystemInfo.WorkingPoint.mainFrequency, 'field', nmrg.SystemInfo.MainField.bz0);
        end
        
        %%
        function obj = dataAcquisition(obj)
            nmrg = obj.boss;
            param = obj.dataPackage.parameters;
            
            try
                freq131 = nmrg.SystemInfo.Xe131FID.f131;
            catch
                error('Cannot get Central frequency.');
            end
            tm = timer('ExecutionMode', 'singleShot', 'StartDelay', param.preDrivingTime);
            tm.TimerFcn = @(myTimerObj, thisEvent) drvXe(nmrg, myTimerObj, freq131 + param.detuning, param.amplitude, param.calibrationField);
            tm.BusyMode = 'queue';
            start(tm);
            
            measurement_data.start = datetime('now');

            nmrg.rbSignal.unsubscribe_signal();
            nmrg.rbSignal.subscribe_signal(ziDemodTime.R, [2 1 1]);
            nmrg.rbSignal.subscribe_signal(ziDemodTime.Theta, [2 1 2]);
            
            measurement_data.T1rhoSignal = nmrg.rbSignal.getDemodSignal(param.duration, 'rate', param.rate);
            
            nmrg.by131.disableDriver();
            measurement_data.ExchangePumpSignal = nmrg.rbSignal.getDemodSignal(param.durationDecay, 'rate', param.rate);
            
            measurement_data.CalibrationSignalTiming = tm.UserData;
            
            measurement_data.finish = datetime('now');
            
            delete(tm);
            obj.dataPackage.measurement_data = measurement_data;
        end
        
        %%
        function obj = analyze(obj)
            nmrg = obj.boss;
            param = obj.dataPackage.parameters;
            measurementData = obj.dataPackage.measurement_data;
            
            signalT1rho  = measurementData.T1rhoSignal;
            singalPump = measurementData.ExchangePumpSignal;

            %---------------------------------------------
            RbIsotopeGyromagneticRatio = 4.7; % Hz/nT for Rb85
            calibrationCoeff = nmrg.SystemInfo.WorkingPoint.mainFrequency/RbIsotopeGyromagneticRatio/nmrg.SystemInfo.MainField.bz0*param.calibrationField;
            [fitting.calibrationValue, fitting.data_t1rho] = findTrigger(signalT1rho, calibrationCoeff, measurementData.CalibrationSignalTiming);
            [fitting.result, fitting.gof, fitting.envelope] = fitenvelopt(fitting.data_t1rho, param.sign);
            [fitting.resultPump, fitting.gofPump] = fitExchangePumping(singalPump, param.sign);

            %---------------------------------------------
            targetResolution = 0.01; % Hz
            L = 2^nextpow2(param.rate/targetResolution);
            fList= (0:L/2)*param.rate/L; Y = fft(fitting.data_t1rho.y-mean(fitting.data_t1rho.y), L);
            P2 = abs(Y/L); P1 = P2(1:L/2+1); P1(2:end-1) = 2*P1(2:end-1);
            
            fitting.fft.fList = fList;
            fitting.fft.value = P1;
            [~, fitting.fft.RabiFreq] = findpeaks(fitting.fft.value, fitting.fft.fList, 'NPeaks', 1, 'SortStr', 'descend');
            
            %-----------------------------------
            obj.dataPackage.analysis_result = fitting;
        end
        
        %%
        function obj = postprocess(obj)
            nmrg = obj.boss;
            param = obj.dataPackage.parameters;
            measuredData = obj.dataPackage.measurement_data;
            analysis = obj.dataPackage.analysis_result;

            %-----------------------------------
            summary.start = datestr(measuredData.start);
            summary.T1rho = analysis.result.tau;
            summary.XeField = analysis.result.a;
            summary.PumpTime = analysis.resultPump.tau;
            summary.finish = datestr(measuredData.finish);
            
            obj.dataPackage.summary = summary;
            nmrg.SystemInfo.Xe131T1Rho = summary;
            
            %-----------------------------------
            getSession;
            sess.addData('GetXe131T1Rho', obj.dataPackage, param.description);
        end
        
        %%
        function obj = presentation(obj, viewer)
            param        = obj.dataPackage.parameters;
            measuredData = obj.dataPackage.measurement_data;
            analysis     = obj.dataPackage.analysis_result;
            summary      = obj.dataPackage.summary;
            
            drvSignalR = measuredData.T1rhoSignal.getNodeData(1);
            drvSignalTheta = measuredData.T1rhoSignal.getNodeData(2);
            pumpSignalTheta = measuredData.ExchangePumpSignal.getNodeData(2);
    
            if nargin == 1
                fig = figure( 'Name', sprintf('Xe131 T1Rho [%s]', datestr(measuredData.start)), 'WindowStyle', 'Docked'); 
                ax1 = subplot(3, 2, 1, 'Parent', fig);
                ax2 = subplot(3, 2, 2, 'Parent', fig);
                ax3 = subplot(3, 2, 3, 'Parent', fig);
                ax4 = subplot(3, 2, 4, 'Parent', fig);
                ax5 = subplot(3, 2, [5 6], 'Parent', fig);
            else
                Table1.NMRGJobs.NMRGJob.updateViewer(viewer, param, summary);
                ax1 = subplot(3, 2, 1, 'Parent', viewer.gui.uiHandle.panelPlot);
                ax2 = subplot(3, 2, 2, 'Parent', viewer.gui.uiHandle.panelPlot);
                ax3 = subplot(3, 2, 3, 'Parent', viewer.gui.uiHandle.panelPlot);
                ax4 = subplot(3, 2, 4, 'Parent', viewer.gui.uiHandle.panelPlot);
                ax5 = subplot(3, 2, [5 6], 'Parent', viewer.gui.uiHandle.panelPlot);
            end
           
            sign = param.sign;
            
            plot(ax1, drvSignalR.x, drvSignalR.y, 'r.-')
            legend(ax1, 'RF Magnetometer demodR', 'Location', 'best');
            xlim(ax1, minmax(drvSignalR.x(:)'));
            xlabel(ax1, 'Time (s)'); ylabel(ax1, 'demodR (V)');
            grid(ax1, 'on');
            
            plot(ax2, drvSignalTheta.x, drvSignalTheta.y, 'b.-');
            xlabel(ax2, 'Time (s)'); ylabel(ax2, 'demodTheta (degree)');
            legend(ax2, 'RF Magnetometer demodTheta', 'Location', 'best');
            title(ax2, sprintf('Calibration Amplitude = %3.2f V', param.calibrationField));
            xlim(ax2, minmax(drvSignalTheta.x(:)')); grid(ax2, 'on');
            
            plot(ax3, analysis.data_t1rho.x, sign.*analysis.data_t1rho.y, 'b.-', analysis.data_t1rho.x, analysis.result(analysis.data_t1rho.x), 'r-'); grid(ax3, 'on');
            xlabel(ax3, 'Time (s)'); ylabel(ax3, 'Xe Field (nT)'); xlim(ax3, minmax(analysis.data_t1rho.x(:)'));
            title(ax3, sprintf('T_{1\\rho} = %3.1f s, b_{max} = %3.1f nT, b_0=%3.1f nT', summary.T1rho, summary.XeField, analysis.result.c));
            
            plot(ax4, analysis.fft.fList, analysis.fft.value, 'r.-'); grid(ax4, 'on');
            xlabel(ax4, 'Frequency (Hz)'); ylabel(ax4, 'FFT'); xlim(ax4, [0, 3.0*analysis.fft.RabiFreq]);
            title(ax4, sprintf('Rabe Freq. = %3.3f Hz', analysis.fft.RabiFreq));
            
            plot(ax5, pumpSignalTheta.x, sign.*pumpSignalTheta.y, 'b.', pumpSignalTheta.x, analysis.resultPump(pumpSignalTheta.x), 'r-'); grid(ax5, 'on');
            xlim(ax5, minmax(pumpSignalTheta.x(:)')); title(ax5, sprintf('Exchange Pumping \\tau = %3.1f s', analysis.resultPump.tau));
        end        
    end
    %%
    methods(Static)
        function res = compare(objList, varargin)
            p = inputParser;
            p.addParameter('isPlot', true, @islogical);
            p.addParameter('viewer', [], @(x)  isa(x, 'Table1.GUI.DataViewer') || isempty(x));
            p.parse(varargin{:});
            
            res = compare@Table1.NMRGJobs.NMRGJob(objList, 'viewer', p.Results.viewer);

            if p.Results.isPlot
                if isempty(p.Results.viewer)
                    fig = figure('Name', sprintf('Xe131 T1Rho'));
                    ax1 = subplot(3, 2, 1, 'Parent', fig);
                    ax2 = subplot(3, 2, 2, 'Parent', fig);
                    ax3 = subplot(3, 2, 3, 'Parent', fig);
                    ax4 = subplot(3, 2, 4, 'Parent', fig);
                    ax5 = subplot(3, 2, [5 6], 'Parent', fig);
                else
                    ax1 = subplot(3, 2, 1, 'Parent', p.Results.viewer.gui.uiHandle.panelPlot);
                    ax2 = subplot(3, 2, 2, 'Parent', p.Results.viewer.gui.uiHandle.panelPlot);
                    ax3 = subplot(3, 2, 3, 'Parent', p.Results.viewer.gui.uiHandle.panelPlot);
                    ax4 = subplot(3, 2, 4, 'Parent', p.Results.viewer.gui.uiHandle.panelPlot);
                    ax5 = subplot(3, 2, [5 6], 'Parent', p.Results.viewer.gui.uiHandle.panelPlot);
                end
                
                drvSignalR = arrayfun(@(d) d.T1rhoSignal.getNodeData(1), res.measuredData, 'UniformOutput', false);
                drvSignalTheta = arrayfun(@(d) d.T1rhoSignal.getNodeData(2), res.measuredData, 'UniformOutput', false);
                pumpSignalTheta = arrayfun(@(d) d.ExchangePumpSignal.getNodeData(2), res.measuredData, 'UniformOutput', false);
                signList = arrayfun(@(p) p.sign, res.param, 'UniformOutput', false);
                t1rhoData = arrayfun(@(d) d.data_t1rho, res.analysis, 'UniformOutput', false);
                t1rhoFFT = arrayfun(@(d) d.fft, res.analysis, 'UniformOutput', false);
                
                cla(ax1); legend(ax1, 'off'); hold(ax1, 'on');
                cellfun(@(drv) plot(ax1, drv.x, drv.y, '.-'), drvSignalR);
                xlabel(ax1, 'Time (s)'); ylabel(ax1, 'demodR (V)');
                hold(ax1, 'off'); grid(ax1, 'on');
                
                cla(ax2); legend(ax2, 'off'); hold(ax2, 'on');
                cellfun(@(drv) plot(ax2, drv.x, drv.y, '.-'), drvSignalTheta);
                xlabel(ax2, 'Time (s)'); ylabel(ax2, 'demodTheta (degree)'); title(ax2, '');
                hold(ax2, 'off'); grid(ax2, 'on');
                
                cla(ax3); legend(ax3, 'off'); hold(ax3, 'on');
                cellfun(@(data, sign) plot(ax3, data.x, sign.*data.y, '.-'), t1rhoData, signList); 
                xlabel(ax3, 'Time (s)'); ylabel(ax3, 'Xe Field (nT)'); title(ax3, 'Rabi Oscillation');
                hold(ax3, 'off');grid(ax3, 'on');
                
                cla(ax4); legend(ax4, 'off'); hold(ax4, 'on');
                cellfun(@(data) plot(ax4, data.fList, data.value, '.-'), t1rhoFFT); 
                xlabel(ax4, 'Frequency (Hz)'); ylabel(ax4, 'FFT'); title(ax4, 'Rabi Frequency');
                hold(ax4, 'off');grid(ax4, 'on');
                
                cla(ax5); legend(ax5, 'off'); hold(ax5, 'on');
                cellfun(@(data, sign) plot(ax5, data.x, sign*data.y, '.-'), pumpSignalTheta, signList);
                xlabel(ax5, 'Time (s)'); title(ax5, 'Exchange Pumping');
                hold(ax5, 'off'); grid(ax5, 'on');
            end
        end
    end
end

function drvXe(obj, timerObj, freq, amp, calibration)
    tic;
    obj.bz0.setDC( obj.SystemInfo.MainField.bz0 + calibration); pause(2.0);
    t1 = toc;
    obj.bz0.setDC( obj.SystemInfo.MainField.bz0); pause(1.0);
    t2 = toc;
    obj.bz0.setDC( obj.SystemInfo.MainField.bz0 - calibration); pause(2.0);
    t3 = toc;
    obj.bz0.setDC( obj.SystemInfo.MainField.bz0); pause(1.0);
    t4 = toc;
    obj.by131.setAC('frequency', freq, 'amplitude', amp, 'range', 10.0); 
    timerObj.UserData = [t1, t2, t3, t4];
end


function [calibrationValue, data_t1rho] = findTrigger(signal, calibrationField, tocList)
    tList0 = signal.getNodeData(2).x;
    data0 = signal.getNodeData(2).y;
    [~, idx] = findpeaks(diff(data0), 'Npeak', 2, 'SortStr', 'descend');
    t0 = tList0(min(idx)+1);
    tList = tList0 - t0;
    
    idxCalibration = tList>0 & tList<tocList(1);
    calibrationValue = mean(data0(idxCalibration));
    
    idxT1rho = tList >= tocList(4);
    data_t1rho.x = tList(idxT1rho) - tocList(4);
    data_t1rho.y = data0(idxT1rho)/calibrationValue*calibrationField;
end

function [result, gof, env] = fitenvelopt(data, sign)
    if nargin == 1
        sign = 1.0;
    end
    tList = data.x;
    env=envelope(sign*data.y, 3000);

    [xData, yData] = prepareCurveData( tList, env );

    % Set up fittype and options.
    ft = fittype( 'a*exp(-x/tau)+c', 'independent', 'x', 'dependent', 'y' );
    opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
    opts.Display = 'Off';
    opts.Lower = [0 0 0];
    opts.StartPoint = [0.5*max(yData) 0.5*max(yData) mean(xData)];

    % Fit model to data.
    [result, gof] = fit( xData, yData, ft, opts );
end

function [result, gof] = fitExchangePumping(signal, sign)
    data = signal.getNodeData(2);
    tList = data.x;
    sigList = sign.*data.y;

    [xData, yData] = prepareCurveData( tList, sigList );

    % Set up fittype and options.
    ft = fittype( 'a*exp(-x/tau)+c', 'independent', 'x', 'dependent', 'y' );
    opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
    opts.Display = 'Off';
    opts.Lower = [0 -Inf 0];
    opts.StartPoint = [max(yData) min(yData) mean(xData)];

    % Fit model to data.
    [result, gof] = fit( xData, yData, ft, opts );
end
